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Abstract 

Hematopoiesis is a complex biological process that leads to the production and reg- 
ulation of blood cells. It is based upon differentiation of stem cells under the action 
of growth factors. A mathematical approach of this process is proposed to carry out 
explanation on some blood diseases, characterized by oscillations in circulating blood 
cells. A system of three differential equations with delay, corresponding to the cell cycle 
duration, is analyzed. The existence of a Hopf bifurcation for a positive steady-state is 
obtained through the study of an exponential polynomial characteristic equation with 
delay-dependent coefficients. Numerical simulations show that long period oscillations 
can be obtained in this model, corresponding to a destabilization of the feedback regu- 
lation between blood cells and growth factors. This stresses the localization of periodic 
hematological diseases in the feedback loop. 

Keywords: delay differential equations, characteristic equation, delay-dependent coefficients, 
stability switch, Hopf bifurcation, cell population models, hematopoiesis, stem cells. 

1 Introduction 

Hematopoiesis is the process by which erythrocytes (red blood cells), leukocytes (white blood 
cells) and thrombocytes (platelets) are produced and regulated. These cells perform a variety 
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of vital functions such as transporting oxygen, repairing lesions, and fighting infections. 
Therefore, the body must carefully regulate their production. For example, there are 3.5 x 
10"'^^ erythrocytes for each kilogram of body weight, so almost 7% of the body mass is red 
blood cells. The turnover rate is about 3 x 10^ erythrocytes/kg of body weight each day, which 
must be carefully regulated by several O2 sensitive receptors and a collection of growth factors 
and hormones. Although understanding of blood production process evolves constantly, the 
main outlines are clear. 

Blood cells, that can be observed in blood vessels, are originated from a pool of hematopoi- 
etic pluripotent stem cells, located in the bone marrow of most of human bones. Hematopoi- 
etic pluripotent stem cells, which are undifferentiated cells with a high self-renewal and differ- 
entiation capacity, give rise to committed stem cells, which form bands of cells called colony 
forming units (CFU). These committed stem cells are specialized in the sense that they can 
only produce one of the three blood cell types: red blood cells, white cells or platelets. Colony 
forming units differentiate in precursor cells, which are not stem cells anymore, because they 
have lost their self-renewal capacity. These cells eventually give birth to mature blood cells 
which enter the bloodstream. 

One can see that the hematopoiesis process is formed by a succession of complex differ- 
entiations from hematopoietic pluripotent stem cells to precursors. These different differen- 
tiations, occurring in the bone marrow, are mainly mediated by growth factors. They are 
proteins acting, in some way, like hormones playing an activator/inhibitor role. Each type 
of blood cell is the result of specific growth factors acting at a specific moment during the 
hematopoiesis process. 

The red blood cells production, for example, called erythropoiesis, is mainly mediated by 
erythropoietin (EPO), a growth factor produced at 90% by the kidneys. Erythropoietin is 
released in the bloodstream due to tissue hypoxia. It stimulates the erythropoiesis in the 
bone marrow, causing an increase in circulating red blood cells, and consequently an increase 
in the tissue p02 levels. Then the release of erythropoietin decreases and a regulation of 
the process is observed: there is a feedback control from the blood to the bone marrow. 
In extreme situations, like bleeding or moving to high altitudes, where needs in oxygen are 
important, erythropoiesis is accelerated. 

White blood cells are produced during leukopoiesis and the main growth factors act- 
ing on their regulation are Granulocyte-CSF (Colony Stimulating Factor), Macrophage- 
CSF, Granulocyte-Macrophage-CSF, and different interleukins (IL-1, IL-2, IL-6, IL-8, etc.). 
Platelets are mainly regulated by thrombopoictin (TPO), which acts similarly to erythropoi- 
etin. 

The hematopoiesis process sometimes exhibits abnormalities in blood cells production, 
causing the so-called dynamical hematological diseases. Most of these diseases seem to be 
due to a destabilization of the pluripotent hematopoietic stem cell compartment caused by 
the action of one or more growth factors. For erythropoiesis, abnormalities in the feedback 
between erythropoietin and the bone marrow production are suspected to cause periodic 
hematological disorders, such as autoimmune hemolytic anemia. Cyclic neutropenia, one 
of the most intensively studied periodic hematological diseases characterized by a fall of 
neutrophils (white blood cells) counts every three weeks, is now known to be due to a desta- 
bilization of the apoptotic (mortality) rate during the proliferating phase of the cell cycle. 

Mathematical models of hematopoiesis have been intensively studied since the end of the 
1970s. To our knowledge, Mackey [Ulin] proposed the first model of hematopoiesis in 1978 
and 1979. This model takes the form of a delay differential equation. Since then it has 
been modified and studied by many authors, including Mackey. The works of Mackey and 
Rudnicki [TH [TS] and Mackey and Rey [TTl [TH [T3] deal with age-maturity structured models 
of hematopoiesis based on the model of Mackey [S]. Recently, Pujo-Menjouet et al. [T^ 



2 



M. Adimy, F. Crauste and S. Ruan 



Hematopoiesis mediated by growth factors 



have studied the model of Mackey [HI [10] and obtained the existence of periodic solutions, 
with long periods comparing to the cell cycle duration, describing phenomena observed during 
some periodic hematological diseases. However, the influence of growth factors has never been 
explicitly incorporated in these models. 

In 1995 and 1998, Belair et al. [3] and Mahaffy et al. [16] considered a mathematical model 
for erythropoiesis. The model is a system of age and maturity structured equations, which 
can be reduced to a system of delay differential equations. They showed that their model 
fitted well with experimental observations in normal erythropoiesis but they stressed some 
difficulties to reproduce pathological behavior observed for periodic hematological diseases. 

In this paper we consider a system of differential equations modelling the evolution of 
hematopoietic stem cells in the bone marrow, of mature blood cells in the bloodstream and 
of the concentration of some growth factors acting on the stem cell population. A delay 
naturally appears in the model, describing the cell cycle duration. This approach is based 
on the early work of Mackey [3 [ID] and the recent work of Belait et al. [3] and Mahaffy 
et al. [H] dealing with erythropoiesis. Our aim is to show that oscillations in such models 
are mainly due to the destabilization of the feedback loop between blood cells and growth 
factors, causing periodic hematological diseases. 

The paper is organized as follows. We first describe the biological background leading to 
the mathematical model. After showing the existence of a positive equilibrium, we analyze 
its local asymptotic stability. This analysis is performed through the study of a characteristic 
equation which takes the form of a third degree exponential polynomial with delay-dependent 
coefficients. Using the approach of Beretta and Kuang [5], we show that the positive steady- 
state can be destabilized through a Hopf bifurcation and stability switches can occur. We 
illustrate our results with numerical simulations and show that long period oscillations can 
be observed in this model, as it can be observed in patients with some periodic hematological 
diseases. 



In the bone marrow, hematopoietic stem cells are divided in two groups: proliferating and 
non-proliferating, or quiescent, cells. The existence of a quiescent phase, also called Go phase, 
in the cell cycle has been proved, for example, by Burns and Tannock [5]. Quiescent cells 
represent the major part of hematopoietic stem cells, that is about 90% of the hematopoietic 
stem cell population. Proliferating cells are cells actually in cycle: they are committed 
to divide during mitosis after, in particularly, having synthesized DNA. Immediately after 
division, proliferating cells enter the Gq phase where they can stay their entire life. 

We denote by Q{t) and P{t) the quiescent and proliferating cell populations at time t, 
respectively. In the proliferating phase, apoptosis, which is a programmed cell death, controls 
the cell population and eliminates deficient cells. We assume that the apoptosis rate, denoted 
by 7, is constant and nonnegative. In the Go phase, cells can disappear, by natural death, 
with a rate 6. They also differentiate in mature blood cells with a rate g{Q)/Q, where the 
function g is assumed to be nonnegative with g(0) = 0, because no cell can become mature 
when there is no hematopoietic stem cell, and continuously differentiable. Moreover, we 
assume that the function Q i— > g{Q)/Q is nondecreasing for Q > 0, which is equivalent to 



2 The Model 



< 5'(0) < 



g{Q) 



< g'iQ) 



for Q > 0. 



(1) 



Q 



It follows, in particularly, that g is nondecreasing and limg^+oo ^(Q) = +oo. 
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Quiescent cells can also be introduced in the proliferating phase during their life in order 
to ensure the population renewal, at a nonconstant rate (3. It is generally accepted that /3 
depends on the total population of non-proliferating cells OUT]. However, the production 
of mature blood cells is also mediated by growth factors through the stem cell population: 
growth factors induce differentiation and maturation of hematopoietic cells via the stem cell 
compartment. Thus, the dependence of (3 on growth factors must be represented. We assume 
that (3 is continuously differentiable. Moreover, in the particular case of erythropoiesis, f3 is 
an increasing function of the erythropoietin concentration, because a release of erythropoietin 
increases the production of red blood cells. Hence, we assume that P is an increasing function 
of the growth factor concentration, with P{Q, 0) = 0, and a nonincreasing function of the Go 
phase population [9l . 

Thus, the equations modelling the differentiation of hematopoietic stem cells in the bone 
marrow are 



dt 

dP 

dt 



-SQ{t) - g{Q{t)) - /3(g(t), E{t))Q{t) 

+2e-T"/3(Q(t - r), E{t - T))Q{t - r), (2) 
-7^(t)+/3(Q(t),i?(t))g(t) 

-e-'^-pm-r),E{t-T))Q{t-T), (3) 



where E{t) is the growth factor concentration at time t. The parameter r > denotes the 
average time needed by a proliferating cell to divide, that is, t is an average cell cycle duration. 
The last term in equation ^ represents the amount of cells coming from the proliferating 
phase at division. They are in fact quiescent cells introduced in the proliferating phase 
one generation earlier. The factor 2 represents the division of each proliferating cell in two 
daughter cells. 

At the end of their development, precursors give birth to mature blood cells, which are 
introduced in the bloodstream. We denote by M{t) the population of circulating mature 
blood cells. These cells only proceed from Gq cells at the rate g{Q)- Mature blood cells are 
degraded, in the bloodstream, at a rate /i > 0. Red blood cells usually live an average of 120 
days, whereas platelets live about one week and white blood cells only few hours. Mature 
blood cell population satisfies the differential equation 

^ = ~MAf(t)+.g(Q(i)). 

The growth factor concentration is governed by a differential equation with an explicit 
negative feedback. This feedback describes the control of the bone marrow production on the 
growth factor production, explained in the previous section. This control acts by the mean of 
circulating blood cells: the more circulating blood cells the less growth factor produced. We 
denote by / the feedback control. The function / depends on the population of circulating 
cells M and is positive, decreasing and continuously differentiable. Then 

'^ = -kE{t) + f{M{t)), 

where /c > is the disappearance rate of the growth factor. In fact, the action of the mature 
blood cell population on the production of growth factor is not immediate: it is slightly 
delayed, but this delay is negligible in front of the cell cycle duration, so we do not mention 
it. 
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At this point, one can notice that system ([l])-® is not coupled: the population in 
the proliferating phase is not needed in the description of the hematopoiesis process, since 
circulating blood cells only come from quiescent cells. From a mathematical point of view, 
the solution of ^ does not depend on the solution of Q whereas the converse is not true. 
Consequently, we concentrate on the system of delay differential equations 

^ = -6Q{t)-9{Q{t))-l3{Q{t),E{t))Q{t) 
+2e-^"/3(Q(t - r), E{t - T))Q{t - r), 
^ ^ -tiM{t)+g{Q{t)), 

^ = -kE{t) + f{M{t)). 

First, one can notice that the solutions Q{t), M(t) and E{t) of system (|4]) are nonnegative. 
Let us suppose, by contradiction, that there exists to > and £ > such that Q{t) > 
for t < to, Q{to) = and Q{t) < for t G {to, to + e). Then, since g{0) = 0, 

^(<o) = 2e-^^l3{Q{to ~ r),E{to ~ r))Q(to - r) > 0. 
at 

This yields a contradiction. Hence Q{t) is nonnegative. Since the functions g and / are 
nonnegative, we similarly obtain that M{t) and E{t) are nonnegative. 

Secondly, one can notice that the solutions of (j4|) are bounded when 6 + g'{0) > 0. 

We first concentrate on the solution E{t). Using a classical variation of constant formula, 
we obtain, for t > 0, 

ft 



E{t) = e-'=*£;(0) + e-"' / e"' f{M{s))ds. 
Jq 

Since the function / is decreasing and bounded, we have 

E{t) < e-''*E{0) + 4^(1 - 6-'=*). 



This yields that, if E{0) < /(0)/fc, then E{t) < f{0)/k and, if /(0)/fc < E{0), then E{t) < 
E{0). Consequently, E{t) is bounded. 

Now we focus on the solution Q{t). If Q is bounded then the mapping t t-^ g{Q{t)) is 
bounded so we will obtain that M{t) is bounded using similar arguments than for the above 
case. 

Let C > be a bound of E and assume that limg^tx; E) — 0, for all -E > 0, and that 
<5 + g'(0) > 0. If 2e-T^/3(0, C) > S + g'{0), then, since the mapping Q P{Q, E) is decreasing 
and tends to zero at infinity, there exists Qq > Q such that 2e^'''^(3{Qo,C) = 6 + g'{0) and 
2e—'^l3{Q, C) <d + g'{0) for Q > Qq. If 2e-'^^/3(0, C) > S + g'{0) then this resuh holds for 
Qo = 0. 

We then set 

Let Q > Qi be fixed and let < y < Q. If y < Qo, then 

2e-^^/3(y,C)y < 2e-^^/3(0, C)Qo - (S + g'mQi < (<5 + g'(0))Q. 
On the other hand, ii y > Qq, then 

2e-^^/3(y, C)y < {6 + g'{0))y < {5 + .g'(0))Q. 
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Thus, 



2e-T" max (3{y, C)y < [5 + g'{0))Q, for Q > Qi. 

0<y<Q 



We assume now, by contradiction, that limsup Q(t) — +00. Then there exists to > t such 
that 

Qit) < g(to), for te[to~ T,to], and Q(io) > Qi- 
Since the function E i— > (3{Q, E) is increasing, we deduce, from JT]), that 

Q'[h) < g'{0)Q{to) - g{Q{to)) - /3(Q(io), S(to))Q(io) < 0. 

We obtain a contradiction so Q is bounded. 

These results are stated in the fohowing proposition. 

Proposition 2.1. Assume that limg^oo f3{Q, E) — 0, for all E > 0, and that S + g'{Q) > 0. 
Then the solutions of system are bounded. 

A solution {Q, M, E) of (j4|) is a steady-state or equilibrium if 

dQ _d'M _dE _^ 
dt dt dt ' 

that is _ _ 

SQ + 9{Q) = {2e-'- -1)(3{Q,E)Q, 

mM = g{Q), (5) 
kE = f(M). 

We make some remarks. It would be nonsense to suppose that the rates /i and k may vanish, 
because we cannot allow the blood cell population to grow indefinitely and growth factor is 
necessarily degraded while in blood. Hence, we assume that ^ > and A; > 0. 

Since 17(0) = 0, it follows that (0, 0, /(0)/fc) is always a steady-state of (|4|) that we will call 
in the following the trivial equilibrium of (|3]) . This steady-state corresponds to the extinction 
of the cell population with a saturation of the growth factor concentration. 

From now on and in the following, we assume that 



,'T,j{Q^-k^[r,^^Q)))-'- 

Since limg^+oo g{Q) = +00 and YvaiM^+oa f{M) = 0, this property holds true if P(Q, 0) = 0, 
for aU Q > 0, or limg^+oo PiQ, E) = 0, for all E>Q. 

Let us assume that dU has a nontrivial positive steady-state {Q* , M* , E*), that is, Q* > 0, 
M* >0 and E* > satisfy ©. Then 

(2e-^^-l)/3(g*,i?*)=<5 + ^^, Ar^^giQ*) and E*^^f{M*). (7) 
Necessarily, we have 

2e"'^"' - 1 > 0, 

which is equivalent to 

. < i^. (8) 

7 
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We assume that dH) holds. Since Q* > and 



= If (AH - If 



we must have 



Let us define 



(2e-^^-l)/3(^Q*,i/Q5(Q*) 



:=/3(^g,^/Q.g(Q))) . (10) 

Since / is decreasing, g is nondecreasing and /3{Q,E) is decreasing with respect to Q and 
increasing with respect to E, we obtain 

'^'<«^i {») - f^/' {» i (^''«)) < " 

Therefore, /3 is decreasing, 

/3(0) = /? ( 0, i/(0) ) and hm /3(Q) = 0. 

Moreover, from ([1]), the function Q i-^ 5 + g{Q)/Q is nondecreasing. Consequently, equation 
([5]) has a positive solution, which is unique, if and only if 



J + .g'(0)< (2e-^--l)/3(^0,i/(0)). 



(11) 



These results are summarized in the next proposition. 
Proposition 2.2. Assume that /i > and k > 0. 
(i) If 

(12) 



^ + g'(0)> (26-^^-1) /3(^0,i/(0)), 
then system 0) /las a unique steady-state (0,0, /(0)/fc); 



(ii) If condition 111]) holds, then system has two steady- states: a trivial one (0, 0,/(0)/fc) 
and a nontrivial positive one {Q*,AI*,E*), where Q* is the unique positive solution of 
{3), M* = g{Q*)/ti and E* = f{M*)/k. 

The above proposition indicates that system Q undergoes a transcritical bifurcation when 
J + g'(0) = {2e^'r^ - 1) /3(0, /(0)/fc). Since the trivial steady-state (0,0, /(0)/fc) corresponds, 
biologically, to the extinction of the cell population and a saturation of the growth factor 
concentration, it is not the more interesting equilibrium. It describes a pathological situation 
that can only lead to death without appropriate treatment. Therefore we focus on the local 
stability analysis of the other steady-state {Q*, M* , E*). One can check that condition (fTTj) . 
which ensures the existence of this steady-state, is equivalent to 

J-}-g'(0)</3 0,^ and 0<r<w:=-ln ^ , ■ (13) 



One can notice that, using the implicit function theorem, we can easily show that the 
steady-states Q*, M* and E* are continuously diffcrcntiable with respect to the cell cycle 
duration r. Moreover, Q* and M* are decreasing functions of t and E* is an increasing 
function of r, with lim^^^„_ (Q*(r), M*(t), (r)) = (0, 0, /(0)/fc). 
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3 Local Stability Analysis 



We concentrate on the study of the stabihty of the nontrivial equihbrium {Q* , M* , E*). 
Hence, we assume throughout this section that > and condition holds. 

We investigate the local asymptotic stability of the steady-state {Q* , M*, E*). The delay 
is often seen as a destabilization parameter (see, for example, Metz and Diekmann JLTJ, Fortin 
and Mackey [6]), so we perform this stability analysis with respect to the delay parameter r, 
which represents the cell cycle duration. We recall that Q* , M* and E* satisfy ([7|). 

We linearize (|H) around the equilibrium {Q* , M* , E*). Set 



q{t) = Q{t) - Q* 
The linearized system is 



m(t) = M{t) - M* and e{t) = E{t) - E* 



dq 
'dt 
dm 
'dt 
de 
'dt 



-Aq{t) + Bq{t - r) - Ce(t) + De{t - r), 
-^m(t) + Gq{t), 
~ke{t) ~ Hm{t), 
where the real coefficients A, B, C, D, G and H are defined by 



(14) 



(15) 



A = S + g'iQ*) + PiQ*,E*) + P[iQ*,E*)Q*, 
B = 2e-'-[P{Q*-E*)+(i[{Q\E*)Q*], 
C = l3!,iQ*,E*)Q* >0, 
D ^ 2e-'^^/3^(Q*,£;*)Q* > 0, 
G = g'{Q*)>0, 
H = -f'{M*) > 0. 

One can notice that these coefficients depend, explicitly or implicitly, on the parameter r 
through the equilibrium values Q* , M* and E*. However, we do not stress this dependence 
when we write the coefficients. Moreover, from the assumptions on /?, g and /, the coefficients 
C, D, G and H are strictly positive. 

In the above definitions, we have used the notations 



Furthermore, one can notice that 
A-B = g'{Q*)^ 

and 



9{Q*) 

Q* 



83 

and p'^{Q^E):^J^iQ,E). 



(2e-T"-l) (3[{Q*,E*)Q* >0 



D-G= {2e-^^ - 1) l3'2iQ*,E*)Q* > 0. 
System ()14p can be written in the matrix form 

dX 



dt 



= AiXit)+A2X{t-T) 
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with 



/ 











-C 




m{t) 




G 







V 


e{t) I 




K 


-H 


-k 





( ^ 





D 


and A2 — 













V 









Consequently, the characteristic equation of ((H]) is given by 

det {XI ~Ai- ^26"^^) = 0, 

which reduces to 

(A + /i)(A + k){\ + A- Be-^^) - GH [C - De'^^) = 0. (16) 



We recall the following result: the trivial solution of system ([T4| . or equivalently the steady- 
state of system Q, is asymptotically stable if all roots of have negative real parts, and 
the stability is lost only if characteristic roots cross the axis from left to right, or right to 
left, that is if purely imaginary roots appear. 

Remark 1. If we linearize system ^ around its trivial steady-state (0, 0, f{0)/k), we obtain 

A = 5 + 5'(0) + /3(0,/(0)/fc)>0, D = 0, 

B = 2e-^-/3(0,/(0)//c)>0, G = g'{0)>0, 

C ^ 0, H = -/'(O) > 0. 

Therefore, the characteristic equation (|16p of the linearized system, around the trivial steady- 
state, becomes 

(A + ^)(A + fc) (A + A - Be-^^) = 0. 
It follows that A = — /i < 0, A — fc < or 

X + A- Be-^'' ^ 0. (17) 



Studying the sign of the real parts of the roots of pT]) , we obtain the following proposition 
which deals with the local asymptotic stability of the trivial steady-state of 

Proposition 3.1. Assume that /i > and k > 0. If condition holds, then the trivial 
steady-state of system ^ is locally asymptotically stable for all t > 0, and if condition 1111]) 
holds, then it is unstable for all t > 0. 

Proof. First notice that, when r = 0,A = i3 — ylsoA>Oif condition (fTTj) holds and A < 
if condition ([T^ holds. 

Let r > be fixed. Setting 1/ = At, the characteristic equation (fT7|) is equivalent to 



{v + AT)e'' -Bt = 0. 
From [S], we know that Re(z^) < if and only if 

At>-1, At-Bt>0, and < C sin(C) - cos(C), 

where C is the unique solution of 

C = -Artan(C), Ce(0,^). 
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Since A > and t > 0, condition At > — 1 is satisfied. 

If we assume that condition (fT^ holds, then A > B so At — Bt > 0. By contradiction, 
suppose that Bt > ^sin(C) — At cos{(). Then, from the definition of (, 

At 

Bt > 



cos(C) ■ 

Since A > B > 0, it foUows that 

1 

1 > 



cos(C) ■ 

Consequently, cos(C) > and ( G (7r/2, tt). We deduce that tan(^) > so 

-ATtan(C) < < C- 

This gives a contradiction. Therefore Bt < (sm{C,) — Atcos(C), and all roots of (fT7|) have 
negative real parts. The trivial steady-state is then locally asymptotically stable for all r > 0. 

Assume now that condition (fTTj) holds. Then A < B and At — Bt < 0. Consequently, 
for all T > 0, (|17p has roots with nonnegative real parts and the trivial steady-state is 
unstable. □ 

The above results indicate that the trivial steady-state of ^ is locally asymptotically 
stable when it is the only equilibrium and unstable as soon as the nontrivial equilibrium 
exists. When the transcritical bifurcation occurs, that is when the two steady-states coincide, 
the trivial steady-state is stable. 

We now return to the analysis of the local asymptotic stability of the nontrivial steady- 
state {Q*,M*,E*) of system Q. 

Equation (fTC|) takes the general form 



with 



P(A,r)-HQ(A,T)e-^" -0 (18) 



P(A,t) = + ai(T)A2 +a2(r)A + a3(T), 
Q{\t) = ai{T)\^ + az{T)\ + aQ{T), 



where 



ai(T) = ^ + k + A, cLi{T) = —B, 

02(7-) = fik + A{fi + k), a^{T) = -B{fi + k), 

asiT) = fikA-GHC, ae{T) = -B^ik + GHD. 
We can check that, for all r G [0, Tmax), 

ai(T) -t- a4(T) =fi + k + A- B>0, 

02(1-) + a5(T) = nk + {A - B){^ + k) > 0, 

and 

asM + aeiT) = fik{A - B) + GH{D - C) > 0. 
We will remember, in the following, that 

a,(T) + a,+3(T) > for i = 1, 2, 3. (19) 
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Let us examine tlie case r = 0. This case is of importance, because it can be necessary 
that the nontrivial positive steady-state of Q is stable when r = to be able to obtain the 
local stability for all nonnegative values of the delay, or to find a critical value which could 
destabihze the steady-state (see Theorem 13. 2p . 

When T = 0, the characteristic equation (fT8|) reduces to 

A''' + [ai(0) + a4(0)]A2 + [02(0) + a5(0)]A + 03(0) + ae{0) = 0. (20) 

The Routh-Hurwitz criterion says that all roots of ([20]) have negative real parts if and only 
if 

ai(0) + a4(0) > 0, 
03(0) + 06(0) >0, 

and 

[ai(0) + a4(0)][a2(0) + 05(0)] > 03(0) + ae{0). (21) 

From (Uni), it follows that all characteristic roots of ([^0]) have negative real parts if and only 
if (HI]) holds. 

Proposition 3.2. When t = 0, the nontrivial steady-state {Q* , M* , E*) of ^ is locally 
asymptotically stable if and only if 

{fi + k) [fik +iA- B){fi + k + A-B)]> GH{D - C), (22) 

where A, B, C , D, G and H are given by il5\) . 

In the following, we investigate the existence of purely imaginary roots A = zw, ciJ G M, 
of (fT5|) . Equation (fT5| takes the form of a third degree exponential polynomial in A. In 
2001, Ruan and Wey [5D] gave sufficient conditions for the existence of zeros for such an 
equation, but only in the case where the coefficients of the polynomial functions P and Q 
do not depend on the delay r, that is when the characteristic equation (|18p is given by 
P(A) -I- Q{X)e~'^'^ — 0. Since all the coefficients of P and Q depend on r, we cannot apply 
their results directly. In 2002, however, Beretta and kuang [4j established a geometrical 
criterion which gives the existence of purely imaginary roots for a characteristic equation 
with delay dependent coefficients. We are going to apply this criterion to equation (fTS]) in 
order to obtain stability results for equation (fT4|) . In the following, we use the same notations 
as in Beretta and Kuang [3] to make things easier for the reader. 

We first have to verify the following properties, for all t e [0, Tmax)' 

(i) P{0,t)+Q{0,t)^0; 

(ii) P{iLO,T) + Q{iuj,T) j^O; 

; |A| ^ oo,ReA > o| < 1; 
(iv) F{llj,t) :— \P{iuj,T)\'^ — \Q{iuj,T)\'^ has a finite number of zeros. 

Properties (i), (ii) and (iii) can be easily verified. Let r e [0,T,nax)- Using ITOl) . a simple 
computation gives 

P(0, t) + Q(0, t) = asir) + aeir) > 0. 

Moreover, 

P{iuj,T) + QiiLU,T) = [-{ai{T) + a4{T))uj'^ + asir) + aeir)] 

+i [-uj^ + (a2(r) + a5(T))tj] , 



(iii) limsup| 



P(A, 
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so (ii) is true. Finally, 



Q(A,t) 






P(A,r) 


A|— >oo 


X 



therefore (iii) is also true. 

Now, let F be defined as in (iv). Since 



and 

we have 
with 



|P(zcj,t)|2 =to^+ [a?(T) - 2a2(T)] lu'^ + [aUr) - 2ai(T)a3(T)] uj^ + aUr) 

\Q{iu, r)p = al{T)u;^ + [aUr] - 2ai{T)ae{T)] lu^ + aUr), 

F{uj,t) = uj'^ + bi{T)uj^ + b2{T)uj^ + bair) 

&i(t) = a2(T) -2a2(T) -a|(T), 

&2(-r) = a|(T) + 2a4(T)a6(r) - 2ai(T)a3(T) - a|(T), 

h{r) = aUr)-aUT). 



One can check that 
and 



&i(t) =^^ + P + A'^ -B^ 



&3(t) = fi^P{A^-B^) + G^H^{C^~D^) + 2fikGH{BD-AC), 

where A, B, C, D, G and H are given by pS]) . It is obvious that property (iv) is satisfied. 

Now assume that A = iw, a; £ R, is a purely imaginary characteristic root of ()18p . 
Separating real and imaginary parts, we can show that (w,t) satisfies 

— ai{T)uj^ + a^{T) = — [— a4(T)a;^ + a6(T)] cos(wr) — a5(r)a; sin(wr), (23) 
—up + a2{T)uj — — a5(T)ijj cos(tJT) + [— a4(T)tj^ + a6(r)] sin(tJT). (24) 



One can check that, if (w, t) is a solution of system then (— w, r) is also a solution 

of ([23|) -(|24 p . Hence, if iuj is a purely imaginary characteristic root of (fT8|) . its conjugate has 
the same property. Consequently, we only look in the following for purely imaginary roots of 
psp with positive imaginary part. 
System (US])- ((Ml) yields 



cos(wt) 
sin(a;r) 



(as — 0104) oj^ + (aiflg + 0304 — 0205) cj^ — 0306 



i|a;"* + {a1 - 20405) + a| 
040;^ + (oifls — 0204 — oe) ti-''^ + (0205 — 0305) a; 
a\uj^ + {a\ — 20406) cj-^ + Og 



(25) 
(26) 



where we deliberately omit the dependence of the o^ on r. 

A necessary condition for this system to have solutions is that the sum of the squares of 
the right hand side terms equals one. By remarking that system can be written 



cos(a;T) = Im 
then this condition is 



(P{ilo,t) 



and sin(a;T) = —Re 



Q{iuj,T)J ' 



12 



M. Adimy, F. Crauste and S. Ruan 



Hematopoiesis mediated by growth factors 



That is 

F{uj,t) = 0. 
The polynomial function F can be written as 

F{uj,T)^h{LO^T), 

where /i is a third degree polynomial, defined by 

h{z, t) := + 6i(t)z2 + 62(r)z + bair). (27) 

We set 

AiT)=bl{T)-3b2iT), (28) 

and, when A(r) > 0, 

^o(t) = ^ . (29) 

We then have the following lemma (see [20]). 

Lemma 3.1. Let r G [0,Trnax) cind A(t) and zq{t) be defined by i28\) and i29\) . respectively. 
Then h{-,T), defined in |^7p , has positive roots if and only if 

bsir) < or b^ir) > 0, A(t) > 0, zo(t) > and h{za{T),T) < 0. (30) 

Proof. Details of the proof are given in [20 , Lemma 2.1. □ 

Conditions A(t) > 0, zo{t) > and ft,(zo(T),r) < cannot be easily checked. In the 
following lemma we express them using the coefficients bi, i — 1,2,3. 

Lemma 3.2. Let t > be such that 63(r) > 0. Then A(t) > 0, zo{t) > and /i(zo(r), r) < 
if and only if 

(i) 62(t)<0 or &i(r) < < 62(t) < ^i^, and 
ill) 2A(t)zo(t) + fei(T)&2(T) - 9&3(r) > 0. 

Proof. Let r be given such that 63 (r) > 0. We do not mention, in the following, the depen- 
dence of the coefficients &i on r. 
We have 

A > if and only if bf > 3&2- 

If 62 < 0, this result holds true. Otherwise, it is necessary that fej > 362- In this latter case, 
if bi < 0, then zq > and, if bi > 0, then zq > if and only if 62 < 0. Therefore zq > if 
and only if 

52 

62 < or &i < < 62 < ^- (31) 
Under the assumption (j3ip . h', given by 

h'{z) ^3z'^ + 2biz + b2, 

has two roots, 

z- = --(5i + d) and z+ = --(&i - d) 
o o 
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wi 



th Z- < z+ and d = \/bi — 3&2 (in fact z+ = zo > 0). A simple computation gives 



Noticing that 

we obtain 

Moreover, 
So 



Mz+) = |(6?-d^)-^ + 63. 



bl~d^ ^ {bi - d){2bl - 362 + bid) = ~3z+{bl + bid + d^), 
2 

h{z+)<0 ^ -z+{bl + bid + d^)+bib2-3b3>0. 

bl + bid + d'^ = d^ + bi{bi +d)=d^ - 3&iz_. 
2 o 

h{z+) < <^ ^+ ^ 25iz+z_ + 6162 - 363 > 0. 

o 



Since zj^z- = &2/3, we eventually obtain 

h{z+)<Q ^ 2rf2z+ +6162 -963 > 0. 

This ends the proof. □ 

From the previous lemma, condition (|30p is equivalent to 

63 (r) < or 63 (r) > and hold true. (32) 

Let us show on an example that condition ([5^ is satisfied. 
Note that 63 can be expressed as 

63(r) = ^l'^k^^{A~B){A + B) + G'^H'\C-D){C + D) 
+2fikGHiB{D ~C) + C{B - A)), 



where A, C, D, G and _ff are defined by ([T5]). Since C - £> < and B - A < 0, then 
63(r) < if ^ + B < and B < 0. Moreover, from the definition of B, it follows that B <0 
ii A + B < 0. Consequently, a sufficient condition for 63 (r) <0\sA + B<0. 

Let us assume that 17 is a linear function, given by g{Q) = GQ, with G > 0, and that 
f3{Q,E) = /3i(Q)/?2(S), with Pi{Q) = 1/(1 + Q"), 71 > 0, and P2 an increasing function 
satisfying /32(0) = 0. Then, for r = 0, 

A + B = [4/3i(Q*) + 3/3UQ*)Q1/32(^*). 
Since E* > 0, /32(-E*) > and A + S < if and only if 

4/?i(Q*) + SmiQ* = ^^(|+"(Q")"y" - °' 



that is 



4 / ^ \iA 

n > - and Q* > 



3 V3n--4 
Let S and G be such that 
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where /? is defined by pH)) . and let n > 4/3 be the unique solution of 

Sn 



Then, for n > n w 6.12, 



3n — 4 V 3n 



l/r, 



In = 



^ Urn) l>^«>^+G' 



so Q* > (4/(3ri - 4))!/" and it follows that A + B <0. 

Consequently, 63(0) < and, using the continuity of 63 with respect to t, we deduce that 
there exists r > such that ([32|) is verified for r € [0, r). 

When the reintroduction rate (3 only depends on the growth factor concentration E, 
condition ([5^ is also satisfied for r close to zero. This is numerically obtained in Section [H 

We set I := [0,r) an interval in which ((5^ is satisfied, with < r < Tmaa;- From the 
above remarks, we can find functions /3, g and /, and parameters values such that r exists. 
For re/ there exists at least uj = lo{t) > such that F{u!{t), t) — 0. 

Then, let ^(r) £ [0, 27r] be defined for t G / by 



cos(6l(r)) 



sin(6l(T)) 



(as - 0104) uj^ + (aiag + 0304 - 0205) - Q3«6 
altj'* + (a| — 20406) + Og 

ajUJ^ + (aios - 0204 - ae) + (Q2Q6 - ^3^5) ^ 
a|w'' + (flg — 20406) + a| 



where w = ^(''')- Since ^(^(t), t) for r G /, it follows that 6 is well and uniquely defined 
for aU T e I. 

One can check that zo;*, with w* = w(t*) > is a purely imaginary characteristic root of 
(fT5| if and only if r* is a root of the function Sn, defined by 

5„(^)=^_^ill±|!^, re/, withneN. 
The following theorem is due to Beretta and Kuang [3]. 

Theorem 3.1. Assume that the function S'„(t) has a positive root r* £ / for some n G N. 
Then a pair of simple purely imaginary roots ±iw(r*) of m8\) exists at t — t* and 



sign • 



dRe{\) 



dr 



sign|^(a;(r ),-)|sign|-^ 



(33) 



Since 



condition ^6'6\ is equivalent to 



. f dRe{\) 
sign ^ 



I dr 



sign 



92 



{u;\r*),T*) \ sign 



dSnir) 



dr 



We can easily observe that S'„(0) < 0. Moreover, for all t & I, 5„(t) > S'„+i(t), with 
n G N. Therefore, if Sq has no root in /, then the Sn functions have no root in / and, if 
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the function 5„(t) has positive roots t E I for some n G N, there exists at least one root 
satisfying 

^(r)>0. 

Using Proposition 13.21 we can conclude the existence of a Hopf bifurcation as stated in the 
next theorem. 

Theorem 3.2. Assume that fi, k > 0, condition ill]) is satisfied and holds true. 

(i) // the function Sq{t) has no positive root in I, then the steady-state (Q* , M* , E*) is 

locally asymptotically stable for all t > 0. 

(ii) // the function So{t) has at least one positive root in I, then there exists t* E I such 

that the steady-state (Q* , M* , E*) is locally asymptotically stable for < r < r* and 
becomes unstable for t > t* , with a Hopf bifurcation occurring when t ~ t* , if and 
only if 



In the next section, we illustrate the results established in Theorem 13.21 We show, in 
particularly, that our model can exhibit long periods oscillations, compared to the delay, 
that can be related to experimental observations in patients with periodic hematological 
diseases. 



4 Numerical Illustrations: Long Periods Oscillations 

Let us assume that the introduction of resting cells in the proliferating phase is only triggered 
by the growth factor concentration E{t), that is (3 = (3{E{t)). This assumption is based on 
hypothesis made by Belair et al. [S] for an erythropoiesis model. It describes, for example, 
the fact that the cell population may only react to external stimuli and cannot be directly 
sensitive to its own size. We assume that fi is given by 

PiE)=(3o^^, /3o>0. 
The functions g and / are defined by 

g{Q) = GQ with G > 0, 

and 

f{M) = - ^— — , a,K>0,r>l. 

This latter function often occurs in enzyme kinetics. It has been used by Mackey [51 [TU] to 
describe the rate of introduction in the proliferating phase and by Belair et al. [3] to define 
the feedback from the blood to the growth factor production. 

With these choices for the functions /3, g and /, our model involves 10 parameters, 
including the delay r. Most of the values of these parameters can be found in the literature. 
This is the case for the stem cells mortality rates in the bone marrow, 5 and 7, which are given 
by Mackey f5] or Pujo-Menjouet and Mackey [H], and for the coefficients of the function / 
and the disappearance rate k, given by Belair et al. [3] and Mahaffy et al. [16] . However, it 
is not so easy to determine values for the other parameters. 
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Experimental data indicate similar values for /i and G, since blood cells are almost pro- 
duced at the same rate than they are lost. Then, we will choose 

^ = 0.02d"^ and G = 0.04 d"^ (34) 

The coefficient /3o represents the maximum rate of introduction in the proliferating phase 
and also the value of /3'(0). It strongly depends on the nature of the growth factor. Using 
data about erythropoiesis, we choose 

f3o = 0.5, (35) 

which is less than the maximal rate of introduction proposed by Mackey [S], but seems 
sufficiently large for erythropoiesis modelling. 
Other parameters are given by (see [31 [S]) 

S = 0.01 d~\ 7 = 0.2 d~\ k = 2.8 d"\ . . 

a = 6570, K = 0.0382 and r^7. ^ ' 

With the above choices for the functions /3, g and /, we can explicitly compute the steady- 
states of system (|4|), Q*, M* and E* . In particularly, one can check that condition ([6]) holds 
true. Condition (fT3ll becomes 



iS + G)ia + k) < (3oa and < r < i In ' ^^"^ 



7 \{S + G){a + k) + Poa 
We set 

a(T) = 2e-^^-l forTG[0,w). 
The bmction a is positive and decreasing on [0, Tmax) and satisfies 

{S + G){a + k) , . ^ , r rn ^ 
'-^ < a(r) < 1 for r e 0, T„,ax)- 

The steady-states of (HJ are then defined by 

_ M 1 faPoa{T)^{S + G){a + k)Y^'' 



G A-iA \^ k{6 + G) 



M* = -Q*, 



/3o«(t) -(S + G)- 
For r e [0,Tmax), Q* and M* are decreasing with 

a^n*<tJ— f aPo-i5 + G)ia + k) \ 
^ ^ ^ - G if k{S + G) J 

and 



K^l^ \ k{5 + G) 
and E* is increasing with 



(3o-{S + G) - k 
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Figure 1: The steady-states Q* (solid line), M* (dashed line), and E* (dotted line) of system 
Q are drawn on the interval [0,Tmax), with T^ax — 2.99, for the parameters values given 
in (|34p . ((55)) and When r — T^ax, E* « 2346 but we have stopped the scale on the 

vertical axis at 30 to improve the illustration clarity. 



For the parameters given in (|34p . ([HS)) and ([M)) . the steady-states are drawn on the interval 
[Q,Tmax) on Figure [T] In this case, 2.99. 

The coefficients A, B, C and D, defined in (fTK)) . become 



A = S + G + (3{E*)>0, C = f3'{E*)Q*>0, 

B = 2e-^^l3{E*) > 0, D = 2e-'^^ (3' {E*)Q* > 0, 

and are all strictly positive. The coefficient G is constant and H is still given by H 
—f'{M*) > 0. One can also check that E* is the unique solution of 

(2e-T^ - l)f3{E*) = S + G. 

Thus, 

A^B=iS + Gf^^^ + ^ 



a(r) 

In particularly, we deduce that 

6i(t) = fi^ + P>0, 

&2(t) = ^i^P + 2GH[G{fi + k + A)-AD], 
&3(t) = GH{D -G){2nkA^GH{G + D)) . 

One can notice that bi is now independent of the delay r. Moreover, since bi > 0, the 
polynomial function h, defined in (I27p . has strictly positive roots if and only if (see Lemma 
Oand Lemma [X^ 63 (r) < or 63 (r) > 0, &2(t) < and 

2A{t)zq{t) + 6i(r)62(r) - 963(r) > 0. 
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Figure 2: CoefScients 62 (t) (left) and 63 (t) (right) are represented for r S [0,T,nax) with 

'^max — 2.99. 



Using Maple 9, we compute the coefficients 62 and 63 for the values in (|35|) and 
Results are presented in Figure [2l Since 63 < on [0,2.92) and 62 is always positive, h has 
positive roots if and only if r e / [0,2.92). In this case, h has exactly one positive root 
for r* e [0,2.92), denoted by z* , and, since /i(0,r) < 0, z* satisfies 



dh 

dz 



(Z*,T*) > 0. 



The function 6*0 is drawn for t E I = [0,2.92) in Figure [31 One can see that there are 
two critical values of the delay r for which stability switches occur. In particularly, from 
Theorem l3.2[ a Hopf bifurcation occurs when r is approximately equal to 1.4. Thus, periodic 
solutions appear. 




Figure 3: Graph of the function S'o(t) for r € [0,Tmaa;) with parameters given by (1^ 
and (|36p . and T^ax ~ 2.99. Two critical values of r, for which stability switches can occur, 
appear. 

In Figure m one can check that Si has no positive root on /. Therefore, there exist only 
two critical values of the delay for which stability switches occur, t = 1.4 and t = 2.82. 
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Figure 4: Graph of the function S'i(r) for the same values as in Figure [3] The function has 
no positive root. 



Using dde23 [53], a Matlab solver for delay differential equation, we can compute the 
solutions of ([H for the above mentioned values of the parameters. Illustrations are showed 
in Figures E] to [51 

Before the Hopf bifurcation occurs, solutions are stable and converge to the equilibrium, 
although they oscillate transiently (see Figure [5]). When the bifurcation occurs, periodic 
solutions appear with periods about 100 days (see Figure [1]). These are very long periods 
compared to the delay r, which is equal to 1.4 day. 

9 r n 2.5 



8- / > 




time t 



Figure 5: Solutions Q{t) (solid line), M{t) (dashed line) and E{t) (dotted hne) of Q are 
asymptotically stable and converge to the steady-state values. Damped oscillations can be 
observed. Parameters values are the same as in Figure [3] with r = 0.5. 

When r increases, periods of the oscillations increase and just before a stability switch, 
for T — 2.82, periods of the oscillations are close to 220 days (see Figure [7]). Then, the 
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Figure 6: When t = 1.4, a Hopf bifurcation occurs and periodic solutions appear, with the 
same period for the three solutions Q{t) (solid Hne), M{t) (dashed line) and E{t) (dotted 
line) of (HJ). Periods are about 100 days. Parameters values are the same as in Figure [31 




Figure 7: For r — 2.8, long period oscillations are observed, with periods close to 220 days. 
Solutions Q{t) (solid line), M{t) (dashed line) and E{t) (dotted line) of Q are unstable. 
Parameters values are the same as in Figure [3] 



5 Periodic Hematological Diseases 

Periodic hematological diseases [7^ represent one kind of diseases affecting blood cells. They 
are characterized by significant oscillations in the number of circulating cells, with periods 
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Figure 8: When r = 2.9, the steady-state is asymptotically stable and solutions Q{t) (solid 
line), M{t) (dashed line) and E{t) (dotted line) of ^ converge to the equilibrium. Parameters 
values are the same as in Figure [3] 



ranging from weeks (19-21 days for cyclical neutropenia [7]) to months (30 to 100 days 
for chronic myelogenous leukemia [7]) and amplitudes varying from normal to low levels or 
normal to high levels, depending on cell types. Because of their dynamic character, periodic 
hematological diseases offer an opportunity to understand some of the regulating processes 
involved in the production of hematopoietic cells. 

Some periodic hematological diseases involve only one type of blood cells, for example, 
red blood cells in periodic autoimmune hemolytic anemia or platelets in cyclical throm- 
bocytopenia [22]. In these cases, periods of the oscillations are usually between two and four 
times the cell cycle duration. However, other periodic hematological diseases, such as cyclical 
neutropenia [7] or chronic myelogenous leukemia [5], show oscillations in all of the circulat- 
ing blood cells, i.e., white blood cells, red blood cells and platelets. These diseases involve 
oscillations with quite long periods (on the order of weeks to months). A destabilization of 
the pluripotential stem cell population induced by growth factors seems to be at the origin 
of these diseases. 

Recently, Pujo-Menjouet et al. [13 [TH] and Adimy et al. PQ[2] considered models for the 
regulation of stem cell dynamics and noticed that long periods oscillations could be observed 
in hematopoiesis models. In this work, we have been able to obtain long period oscillations 
(about 100 days) for a blood production model mediated by growth factors, as a result of 
the feedback from blood to growth factors. This result stresses the role of growth factors in 
the appearance of periodic solutions in hematopoiesis models. 

One can notice that by assuming that the cell cycle duration is constant, values of r for 
which a nontrivial steady-state exists are limited and cannot be too large. This does not 
appear in a model with distributed delay, as studied by Adimy et al. [H [2] • 

Numerical simulations demonstrated that long period oscillations in the circulating cells 
are possible in our model even with short cell cycle durations. Thus, we are able to char- 
acterize some hematological diseases, especially those exhibit a periodic behavior of all the 
circulating blood cells. 
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